Analytic Continuation of Quantum Monte Carlo Data 
by Stochastic Analytical Inference 
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We present an algorithm for the analytic continuation of imaginary-time quantum Monte Carlo 
data which is strictly based on principles of Bayesian statistical inference. Within this framework 
we are able to obtain an explicit expression for the calculation of a weighted average over possible 
energy spectra, which can be evaluated by standard Monte Carlo simulations, yielding as by-product 
also the distribution function as function of the regularization parameter. Our algorithm thus avoids 
the usual ad-hoc assumptions introduced in similar algortihms to fix the regularization parameter. 
We apply the algorithm to imaginary-time quantum Monte Carlo data and compare the resulting 
energy spectra with those from a standard maximum entropy calculation. 

PACS numbers: Insert valid PACS here! 



I. INTRODUCTION 

Quantum Monte Carlo simulations are a powerful 
computational tool to calculate properties of interacting 
quantum manyparticle systems, such as spin models or 
strongly correlated electron systems. Of particular inter- 
est in those systems are dynamical correlation functions 
like single-particle spectra or susceptibilies respectively 
dynamical structure factors. However, QMC presently 
provides data only on the imaginary time axis, and the 
necessary analytic continuation of these data has proven 
to be difficult. 

The standard tool to solve this problem is the Maxi- 
mum Entropy Method (MEM) [T]. It uses arguments of 
Bayesian logic [H |3] to obtain the most probable energy 
spectrum. In order to solve this optimization problem ef- 
ficiently, the maximum entropy method approximates all 
occurring probability distributions to be of a Gaussian 
shape. 

In the past efforts were made to provide an alternative 
to this approach [H EJ E] . It was proposed to perform a 
Monte Carlo average over a wide range of spectra instead 
of selecting a single spectrum. So far, the method lacked 
a rigorous rule to eliminate a regularization parameter 
inherent in the algorithm. Although this approach has 
been interpreted in terms of Bayesian inference , non 
of the authors have utililized Bayesian logic to eliminate 
the regularization parameter. 

We show that this stochastic approach can also be un- 
derstood in terms of Baysian statistical inference. We 
derive a strict criterion to eliminate the free parame- 
ter, that is completely based on Bayesian logic. It uses 
Monte Carlo techniques to both calculate the average 



spectrum and to eliminate the regularization parameter. 
It treats all probabilities exactly and hereby avoids the 
approximations made in the maximum entropy method. 
We apply the algorithm to imaginary-frequency quantum 
Monte Carlo data and compare the resulting spectra with 
results from maximum entropy calculations. 



II. THE PROBLEM OF ANALYTIC 
CONTINUATION 

For a finite temperature T quantum Monte Carlo 
simulations can provide accurate estimates Gn for ei- 
ther imaginary-time correlation function G{t) at a finite 
set of N imaginary-time points t„ or, alternatively, for 
imaginary-frequency correlation functions G(iu)n) at a fi- 
nite set of N Matsubara frequencies w„. The frequencies 
are defined as w„ = (2n + Vjir/ (3 for fermions and as 
ujn = 2711^/(3 for bosons with (5 = 1/k^T. 

Because of the stochastical nature of Monte Carlo al- 
gorithms each of the Gn possesses a known statistical 
error. Moreover, the data for the different time or fre- 
quency points are usually highly correlated. Therefore 
the input to the analytic continuation procedure consists 
of the Monte Carlo estimates Gi and their covariance 
matrix 



Cn 



GnGjn GriGfi 



(1) 



In principle the spectral function A{uj) = — ^ImG(aj 
iO+) can be extracted from these data by inverting 



dwif„(w)A(w) 



(2) 



with 
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1 ± e-'^f^ 



(3) 
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for time dependent data or 
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(4) 



for frequency dependent data, where the upper sign holds 
for ferniions and the lower one for bosons. The spectral 
function is normalized to 



M= dwA{uj) 



(5) 



and is nonnegative for all lj. However, a direct inversion 
of Eq. [2] is an ill-posed problem and numerically impos- 
sible. 

A least-square fit of A{lo) to the data Gn minimizes 
the -estimate 

(6) 

with respect ro ^(1^). This approach leads to a multitude 
of different solutions and consequently cannot solve the 
problem either. 



A. The Maximum Entropy Method 

The maximum entropy method can be understood as 
an attempt to regularize the least-square fit described 
above. One defines the entropy 



S[A] 



dbj A(lo) In 



AM 

D{u;) 



(7) 



relative to a default model D(uj). Any information, that 
is known about the spectrum beforehand, can be encoded 
in the default model. If D{uj) is nonnegative and pos- 
sesses the same norm A/" as the spectrum A{uj), the en- 
tropy 5* will be nonpositiv and maximal for D(uj). In- 
stead of just minimizing the MEM minimizes the 
quantity 



Q[A]^x'[A]-aS[A] 



(8) 



introducing a regularization parameter a. This optimiza- 
tion problem can be numerically solved for fixed a to find 
the minimizing spectrum Aa(uj). In the limit of a — > cxd 
the spectrum minimizing Q is the default model D{lu). 
For a ^ the least-square fit is regained. Thus the 
parameter a interpolates between the fit result and the 
default model. 

In order to find a criterion to eliminate the parameter 
different approaches exist. The simplest rule is to take 
the spectrum where ~ 1. This choice ensures that the 
differences between model and data are of the order of 
the error bars thereby avoiding overfitting. In order to 
derive more sophisticated methods the MEM needs to be 
reinterpreted by means of Bayesian statistical inference 

mm- 



The MEM can be reformulated by defining subjective 
probabilities for the quantities involved in the analytic 
continuation problem. Let P[A] denote the prior proba- 
bility of the spectrum A{uj). P[A\G] denotes the poste- 
rior probability of A given the input data G and ^[6*1 A] 
the likelihood function. Bayes 's Theorem [7j relates these 
probabilities to each other: 



P[A\G] = P[G\A] P[A] I P\G\. 



(9) 



The probability P\G\ is called the evidence and serves as 
normalization for the posterior probability P[A|G']: 



P{G\ = \VAP{G\A\ 



One indentifies 



P\G\A\ = -^M-x\A\) 
^1 



and 



The quantities 



and 



P[A]- — exp(a5[A]). 
^2 



(10) 



(11) 



(12) 



(13) 



(14) 



normalize the respective probabilities. This way the pos- 
terior probability can be rewritten as 



P[A\G] = 



.-QIA] 



with 



Z,Z2P[GY 
Z\Zi 



(15) 



(16) 



Thus the minimization of Q can be reinterpreted as the 
maximization of the posterior probability P[A|G'] ~ e"'^. 
The MEM therefore determines the most probable spec- 
trum AoL given the input data G. 



2. Bayesian Inference and the Regularization Parameter a 

This alternative formulation of the problem provides 
the necessary tools to eliminate the free parameter a [5J 
IH]. Eq. |9]can be rewritten including a: 

P[A,a\G] = P[G\A,a] P[A,a] / P[G]. (17) 



3 



If one applies Bayes's theorem to factorize P[A, a] and 
integrates over A, the relation 



P[a\G] ^P[a] J VAP[G\A,a] P[A\a] / P[G] 



P[a] 



ZiZ2P[G] 



(18) 



for the posterior probabilty P[a|G'] can be found. Anal- 
ogous to the argument given above, one identifies 
P[G\A,a] - exp(-x^[A]) and P[A\a] - exp(aS'[A]). 
The evidence 



P[G\ = / da 



P[a] jVAe 



-Q\A\ 



(19) 



is an Qf-independent normalization constant. All quan- 
tities in this equation are known except P[a], the prior 
probability of a. It is either taken to be constant or to 
be the Jeffreys prior 1/a [HI [101 E] • However, the choice 
of P[a] turns out to be of little influence on the resulting 
spectra. 

By assuming all probabilities involved to be of a Gaus- 
sian shape a numerical treatment of the equations ( 15 ) 
and (18 1 is possible. There are two alternatives: 



1. One calculates a* as the a that maximizes P[a|G] 
and takes A^* as the final result for the spectral 
function 019]. 

2. One averages over all A^ weighted by the posterior 
probability of a, i.e. the average spectrum 



doiP{a\G\Ao, 



(20) 



is taken as the final result [TT] . 



It is not a priori clear, which of the two algorithms is 
favorable. 



B. Stochastic Analytical Inference 

Stochastic Analytical Inference is an alternative to the 
standard MEM which does not employ the explicit regu- 
larization of the fit by the entropy Eq. ([t]) . Rather than 
maximizing P[^|G] an average over all possible spectra 
weighted by 



w ~ exp(— x^/a) 



(21) 



is performed. Beach refined this approach, by introduc- 
ing the default model D{ijj) of the MEM into the algo- 
rithm [5 . By mapping uj unto x G [0, 1] using 



(22) 



a dimensionless field n(x) can be defined: 
_ A{<\>-\x)) 

The field n(x) is normalized to 1: 



(23) 



1 = J dx n{x) . 



By calculating the average field 

{n{x))a = ^ / T^'n{x) n{x) e-x'["(")I/^ 



with 



The measure 



Z= /X''n(a;)e->^'["(^)1/" 



(24) 



(25) 



(26) 



V'n{x) ^Vn{x)Q[n] S\ Jdxn{x)~l\ (27) 

restricts the integration to fields n{x) that satisfy the 
norm rule Eq. (24) and the positivity requirement. In 



Eq. (27) 



1, if Va; : n{x) > 
otherwise 



The average spectrum (A) a can be regained via 
{A{uj))„ = D{Lo){n{cl>{u;)))0 . 



(28) 



If is interpreted as a Hamiltonian of a fictitious physi- 



cal system, Eq. (25) possesses the structure of an canon- 
ical ensemble average at a temperature a. The laws of 
statistical mechanics then state, that the average spectral 
function (A) a minimizes the free energy 



aS 



(29) 



This expression displays a similar structure as Eq. (|8|. 
Thus the averaging process implicitely generates an en- 
tropy S. However, this entropy does not have the explicit 
form of Eq. (ItI . In the limit a — > the averaging process 
minimizes y . Whereas in the limit a ^ oo the aver- 
age in Eq. ( 25 1 is completely uneffected by and will 

(n{x)) 



constrainded by Eq. 24 



result in {n{x)) = 1. In this 
case the resulting spectrum is the default model. The al- 
gorithm therefore exhibits the same limiting cases as the 
MEM. Additionaly, Beach has shown that a mean field 
treatment of the ficticious physical system described by 
X^ is formaly equivalent to the MEM [5|. 

The remaining open question, namely how to eliminate 
the parameter a, was addressed by all preceding authors 
differently. 
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1 . Sandvik proposes to examine the plot of the average 
entropy against a and identifies the final a by a 
sharp drop in the entropy curve [4;. 

2. Beach examines a double-logarithmic plot of the 
average ^^id identifies the the final a by a kink 
in the x^-curve [5]. 

3. Syljuasen argues to take a = 2Nb/N, where Nb is 
the number of Monte Carlo bins and N the number 
of input data points, as before [S]. 

All criteria are merely based on heuristic arguments. The 
simple rule to take ~ 1 is also applicable to this 
method and should be mentioned here. 



1. Bayesian Statistical Inference 

In the following we will use Bayesian inference to de- 
rive a new criterion to eliminate the regularization pa- 
rameter a. In constrast to the MEM the Stochastic Ana- 
lytic Continuation does not maximize the posterior prob- 
ability P[A|G']. Instead, It averages all possible fields n 
(omitting the argument x in the progress) weighted by 
P[n\G]: 



VnnP[,i\G] 



(30) 



Bayes's theorem can be applied to factorize P[7i|G] anal- 
ogous to Eq. (|9|): 



P[n\G\ = P[G\n] P[n] / P[G] 



(31) 



The Stochastic Analytic Continuation does not introduce 
an explicit entropy term. Following Ref. 6 only the pos- 
itivity requirement and the norm rule Eq. ( 24 ) enter the 
prior probability 



P[n] ^ e {n{x)) S \ J dx n{x) - 1 

^0 

The likelihood function is identified as 

1 



P[G\n] = -e 



-X-'/c 



(32) 



(33) 



the posterior probability results in 



P[n\G] = e {n{x)) 5\ I dx n{x) - 1 | -g-x'!"!/" 



(36) 



as expected from the comparison of the equations (25 1 



and (30) 



2. Bayesian Inference and the Regularization Parameter a 

Bayesian logic can also be utilized to calculate the 
posterior probability P[a|G]. Substituting n for A in 
Eq. (18 1 and identifying P[n|Q;] = P[n] with P[n] from 
Eq. (32) and correspondingly P[G|n, a] = P[G'|n] with 
P[G|n]from Eq. (33), one obtains 



P[a|G] =P[a] VnP[G\n,a] P[n\a] / P[G] 



'Z'P[G] 



The evidence 



P[G] = / da 



P[a]e 



Z'{a) 



(37) 



(38) 



is again an a-independent normalization constant. The 
combination of the equations ( 34 1 and ( 37 1 gives the final 



expression for the a-dependence of the posterior proba- 
bility: 



P[a|G] ^ P[a] a 

N/2 



-X [n]/a 



(39) 



Analogous to the MEM one has two possibilities to treat 
the regularization parameter: 

1. One calculates a* as the a that maximizes P[q;|G] 
and takes (n)cc as the final result. 

2. One averages over all {n)a weighted by the poste- 
rior probability of a, i.e. the average field 



{{n)) = J daP[a\G]{n)^ 
is taken as the final result. 



(40) 



By evaluating a Gaussian integral the normalization Z' 
can be readily calculated to be 

Z' = J 2?Ge-x'/" = {2TTa)^^'^VdetC . (34) 



Usins 



P[G] = / P' 



-X^[n]/a ^ 



(35) 



III. MONTE CARLO EVALUATION 
A. Configuration and Update Scheme 

In order to calculate the quantities appearing in equa- 
tions ( 25 1 and ( 39 ) a numerically treatable approxima- 



tion for the field configuration n{x) and the intergration 
measure Vn has to be found. Our implementation closely 
follows Ref. 5 , The field configuration is represented by a 
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superposition of delta function walkers with residues r„ 
and coordinates a;„: 



n{x) = S^TnSix - Xn)- 



(41) 



The Monte Carlo updates consist of randomly proposed 
shifts of the coordinates x„ and random redistributions of 
the residues r„. Redistributions that are not only norm 
conserving but also conserve higher moments of the con- 
figuration [5^ , have proven to be effective as well. 



The average Eq. ( 25 ) is evaluated by a standard Monte 



Carlo simulation using Metropolis weights. The regular- 
ization parameter a is treated as the temperature of the 
system. The simulation is performed for a wide range 
of different a- values. A parallel tempering [T^ [2] 
algorithm is necessary to ensure convergence for small 
a. In order to measure the average field configuration a 
histogram of the delta function walkers is recorded. 



B. Calculation of the probability P[a|G] 

A particular problem in the proposed approach is that 
a numerical treatment of Eq. ( 39 ) involves the calculation 
of the quantity 



Z 



(42) 



This is equivalent to calculating a partition function in a 
canonical ensemble at temperature a. Standard Monte 
Carlo techniques are only able to calculate thermal ex- 
pectation values but not the partition function itself. We 
use a Wang-Landau algorithm [T3 [H] to generate the 
density of states p{E) of the system. Once p{E) is cal- 
culated, the partition function can be obtained by 



Z ^ JdEp{E)e-^/". 



(43) 



The Wang-Landau algorithm performs a random walk in 
energy space with probability p{E) = 1/ p{E) using the 
usual metropolis weights. Since the density of states is 
unknown at the beginning of the simulation, one starts 
with an arbitrary starting value, e.g. p{E) = 1. For 
each visited energy one updates an energy histogram and 
multiplies the density of states by an modication factor 
/ > 1. When the histogram is reasonably flat, one re- 
sets the histogram and restarts the simulation with a new 
modification factor /' — \f] . The starting value of / is 
usually taken to be Euler's constant and the procedure is 
repeated until / is very close to 1 (16 times in our imple- 
mentation). The resulting piE) is the density of states 
of the system up to an unknown normalization factor. 
In order to speed up the convergence of the algorithm, 
it is advisable to divide the energy range of interest into 
several slightly overlapping smaller inter vails. 



IV. SIMULATION RESULTS 

A. The model 

We apply the algorithm to imaginary-time data from 
quantum Monte Carlo simulations. As test case we con- 
sider the two-dimensional single-band Hubbard model 



(44) 



(» .j)o" 



Here i and j are lattice site indices, the operators cj^ {cia) 
create (destroy) an electron with spin a e at site 

i, riia = c\^Ci„ is their corresponding number density, 
t is the hopping parameter between neighbouring sites 
(denoted by (i, j)) and U implements the local Coulomb 
repulsion. The full lattice model was approximated by a 
two by two cluster embedded in a mean field using the 
Dynamical Cluster Approximation [TTJ UHl [19] . Using a 
weak-coupling expansion in continuous imaginary time 
[2ni m] the single-particle Green function 



p 



(45) 



was calculated for a certain number of Matsubara fre- 
quencies ujn = (2n -|- l)7r//3. Here T is the imaginary- 
time ordering operator, (•) denotes a thermal expecta- 
tion value and Ci(T) = e~^'^ CiC^'^ . The model was simu- 
lated for U = W , where W = denotes the bandwidth, 
and a fixed filling (n^) = 0.9 for several temperatures 
T. Within the weak coupling expansion it is possible to 
calculate the Green function directly in frequency space 
|21j . so that no Fourier transformation or discretization 
of the imaginary time axis is necessary. In all simulations 
the number of measured matsubara frequencies was re- 
stricted to rimax — '^U [3, which has proven to be sufficient 
for all calculation. A further increase of the number of 
frequencies had no influence on the analytic continuation 
results. 



B. Monte Carlo results 

Fig [T] shows the a-dependence of the single particle 
spectra calculated by the parallel tempering Monte Carlo 
simulation (/? — 1411^"^). A Gaussian default model 



D{u) 



72cr 



(46) 



with cr = 1 was used. The shape of the default model is 
clearly visible for large a. One can see how several differ- 
ent peaks and other structures appear for for decreasing 
a. Since the a-dependence is so strong one definitely 
needs a criterion to eliminate the regularization parame- 
ter. 
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FIG. 1: Simulated spectra for a range of regularization pa- 
rameters a {P = 1AW~^). For large a the Gaussian shape of 
the default model is visible. For decreasing a several features 
begin to appear. 



The density of states calculated by the Wang-Landau 
simulation and the probability distribution ^[010] fol- 
lowing Eq (39 1 is shown in Fig. |2| P[a|G] is plotted for 



the two most common choices for P[a] , i. e. P[q\ — const, 
and P[a] — 1/a. The density of states varies over at least 
15 orders of magnitude (note the logarithmic scales). The 
probability distributions P[q!|G'] exhibit a well defined 
peak at d ~ 0.03. Note that the two different choices 
for P\a\ have only weak influence on the position of the 
peak. The two different probability distributions are used 
to calculate the final single particle spectrum. Following 
the discussion in section [ll B 2[ Fig. [2] shows the average 
of all spectra of Fig [l] weighted by P[a|G'] and the spec- 
trum whose a maximizes P[a|G']. The resulting spectra 
are nearly indistinguishable and show that neither the 
ambiguity in the treatment of the probability distribu- 
tion nor the choice of P[a\ have a significant infiuence on 
the resulting spectrum. 

Let us compare our results from the stochastic analyt- 
ical inference with those obtained with other methods to 
fix a. Fig. [3] shows that the point where ~ 1 corre- 
sponds to a ~ 0.1 and that the x^-^stimate exhibits a 
kink in the same a-region. The chosen a = 0.1 is larger 
than a ~ 0.03. That indicates that the spectra deter- 
mined with this criterion are stronger regularized than 
the spectra calculated by the probability distributions 
in Fig. [2j However, at least for the QMC data under 
consideration, the difference between the two spectra are 
only small. The spectrum for a — is also shown 
in Fig [3] For the present data set this corresponds to 
a ~ 140 which is far more regularized than the spectra 
determined by all other methods. The entropy (Fig. ^p) 
shows no significant features and gives no indication how 
to choose the a-parameter. A sharp drop in the entropy 
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FIG. 2: The probability distributions P[q]G] (a) based on a 
Wang-Landau simulation of the density of states (b). The 
different choices for P[a] only have a weak influence on the 
position of the peak. The resulting spectra (c) are calculated 
by either averaging all spectra over P[a|G] or by taking the 
spectrum that maximizes them. The four different spectra 
are practically identical. 



curve is not visible in the simulated area. 

Finally, we compare the SAI with the standard MEM 
approach. Fig. [4] shows results of a Maximum Entropy 
calculation using Bryan's algorithm [11] for the same 
QMC data as in the previous section. The qualitative 
behaviour is similar to the SAI simulation: The proba- 
bility distribution P[a|G] shows a noticeable dependence 
on the prior probability P[a]. However, the resulting 
spectra depend neither on P[a] nor on wether one aver- 
ages over P[a|G'] or wether one takes the maximum. The 
X^ ~ 1 rule determines an a which is again larger than 
the one calculated by Bayesian inference. Accordingly 
the spectrum calculated by this criterion is more regu- 
larized, although here the difference is relatively small. 
Interestingly, in MEM the interesting values for a are 
about one or two orders of magnitude larger compared 
to those appearing in the SAI simulations. There seems 
to be no direct correspondence between the a-values of 
the two methods. 

An extended comparison of SAI spectra with results 
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FIG. 3: The double logarithmic plot of (a) shows a kink 
at a ~ 0.1 which is also the region where ~ 1- The 
resulting spectrum (c) is similar to the one in[2jD. The entropy 
(b) shows no significant features and gives at least for this 
data set no indication how to determine a. The spectrum for 
a = 2Nb /N ~ 140 (c) is much more regularized than the 
spectra determined by all other methods. 
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FIG. 4: Results of a Maximum Entropy calculation for the 
same QMC data as in Fig. [2] and Fig. [3] The probability dis- 
tribution P[q|G] (a) shows a noticeable dependence on P[a], 
but analogous to SAI, the resulting spectra (c) are identical. 
The a where x^ ~ 1 (b) is larger (a ~ 100) than the one for 
which P[q|G] is maximal {a ~ 10). Accordingly the spectrum 
chosen by the x^ ~ 1 rule is more regularized than the one 
calculated by Bayesian inference. 



CONCLUSION 



of maximum entropy calculations for several tempera- 
tures is collected in Fig. [5] All calculations are based on 
the Gaussian default model Eq (46). As already noted 
before, MEM tends to stronger regularize the spectra 
and consequently do the SAI spectra exhibit noticeably 
sharper features for all temperatures shown. Especially 
the pseudo-gap, that opens at beta = 3AW~^, is cap- 
tured nicely by SAI while the MEM cannot resolve it yet 
at that temperature. 

One last, but nevertheless important, questions con- 
cerns the dependence of the spectra on the default model. 
To this end we show in Fig.[6]again SAI and the MEM re- 
sults for the spectrum at /3 = 34 W"^, this time however 
based on a different default model, namely a rectangu- 
lar default model of width 3.6. The resulting spectra are 
very similax to the one obtained for the Gaussian default 
model presented in Fig. |5] Thus, even at low tempera- 
tures the resulting spectra are quite independent of the 
default model. More precisely, we could not detect a sig- 
nificant default model dependence at any temperature. 



We have demonstrated that the stochastic analytic 
continuation method introduced by Sandvik and Beach 
can be interpreteted in terms of Bayesian probability the- 
ory. We developed an algorithm that uses Monte Carlo 
techniques to both calculate the average spectrum and 
to eliminate the regularization parameter. It treats all 
probabilities exactly and hereby avoids the approxima- 
tions made in the maximum entropy method. 

Comparisons to different approaches to fix the regu- 
lariazation parameter a and standard MEM show that 
the SAI results in robust spectral functions which are 
less regularized and consequently show more pronounced 
features, in particular with decreasing temperature in the 
model calculations. As known from standard MEM, no 
significant dependence on the default model could be ob- 
served. 

On apparent drawback of the method is the necessity 
to perform simulations for a broad range of values for a. 
Although this can be performed with parallel tempering 
techniques, the required computer resources are orders 
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FIG. 5: Spectra simulated by Stochastic Analytic Inference compared to Maximum Entropy Calculations. All calculations are 
based on a Gaussian default model 
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FIG. 6: A Stochastic Analytic Inference result for /3 = 34 
based on a flat default model. The spectrum is very similar 
to the one shown in Fig. [2j;. We conclude, that for the QMC 
data under consideration the calculated spectra are quite in- 
dependent of the default model. 



of magnitude larger than for standard MEM approaches. 
As the resulting spectra tend to be less regularized one 
has to ponder the gain in details in the structures against 
the dramatic increase in computer time. 
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